12. Self-calibrated CS-pMR image reconstruction from undersampled non-Cartesian data

In this tutorial we will reconstruct an MRI image from radial undersampled kspace measurements. Let us denote $\Omega$ the undersampling mask, the under-sampled Fourier transform now reads $F_{\Omega}$.

We use the toy datasets available in pysap, more specifically a 2D brain slice and under-sampled Cartesian acquisition over 32 channels. We compare zero-order image reconstruction with self-calibrated multi-coil Compressed sensing reconstructions (analysis vs synthesis formulation) using the FISTA algorithm for the synthesis formulation and the Condat-Vu algorithm for the analysis formulation. The multicoil data $(y_\ell)_\ell$ is collected across multiple, say $L$, channels. The sensitivity maps $(S_\ell)_\ell$ are automically calibrated from the central portion of k-space (e.g. 5%) for all channels $\ell=1, \ldots, L$.

We remind that the synthesis formulation of the non-Cartesian CS-PMRI problem reads (minimization in the sparsifying domain):

$$ \widehat{z} = \text{arg}\,\min_{z\in C^n_\Psi} \frac{1}{2} \sum_{\ell=1}^L\|y_\ell - F_\Omega S_\ell \Psi^*z \|_2^2 + \lambda \|z\|_1 $$

and the image solution is given by $\widehat{x} = \Psi^*\widehat{z}$. For an orthonormal wavelet transform, we have $n_\Psi=n$ while for a frame we may have $n_\Psi > n$.

while the analysis formulation consists in minimizing the following cost function (min. in the image domain):

$$ \widehat{x} = \text{arg}\,\min_{x\in C^n} \frac{1}{2} \sum_{\ell=1}^L \|y_\ell - F_\Omega S_\ell x\|_2^2 + \lambda \|\Psi x\|_1 \,. $$

• Author: Chaithya G R & Philippe Ciuciu
• Date: 01/06/2021
• Target: ATSI MSc students, Paris-Saclay University

# 导入 MRI 重建运算符
# Package import
from mri.operators import NonCartesianFFT, WaveletN, WaveletUD2
'''
NonCartesianFFT : 非笛卡尔傅立叶变换（Non-Cartesian FFT），用于处理 螺旋或径向 MRI 采样 的傅立叶变换。
WaveletN : 正交小波变换（如 Symmlet-8 或 Daubechies），用于 合成公式（Synthesis Formulation） 压缩感知重建。
WaveletUD2 : 未降采样双正交小波变换（Undecimated Bi-Orthogonal Wavelet Transform），用于 分析公式（Analysis Formulation） 压缩感知重建。
'''
# 导入 MRI 处理的实用工具
from mri.operators.utils import convert_locations_to_mask, \
    gridded_inverse_fourier_transform_nd
'''
convert_locations_to_mask : 将 k-space 采样坐标转换为二值掩码，用于 重建欠采样的 MRI 数据。
gridded_inverse_fourier_transform_nd : 基于网格插值的逆傅立叶变换，通常用作 MRI 重建的初始估计。
'''
#  导入 MRI 图像重建器
from mri.reconstructors import SelfCalibrationReconstructor
'''
SingleChannelReconstructor : 单通道 MRI 迭代重建器，支持 FISTA、POGM 和 Condat-Vu 等压缩感知优化算法。
'''
from mri.reconstructors.utils.extract_sensitivity_maps import get_Smaps

# 导入 MRI 样本数据
import pysap
from pysap.data import get_sample_data

# Third party import
from modopt.math.metrics import ssim
from modopt.opt.linear import Identity
from modopt.opt.proximity import SparseThreshold
'''
ssim : 结构相似性指数（SSIM），用于 评估重建图像的质量，与原始 MRI 进行比较。
Identity : 单位变换（Identity Transformation），即 不对数据做任何变换，在优化过程中常用于保持数据结构。
SparseThreshold : 稀疏阈值运算（L1 正则化 / 软阈值收缩），用于 压缩感知重建中的稀疏性约束。
'''

import numpy as np
import matplotlib.pyplot as plt

Loading input data

# Loading input data
cartesian_ref_image = get_sample_data('2d-pmri')
#image = pysap.Image(data=np.sqrt(np.sum(cartesian_ref_image.data**2, axis=0)))
image = np.linalg.norm(cartesian_ref_image, axis=0)

# Obtain MRI non-Cartesian radial mask
radial_mask = get_sample_data("mri-radial-samples")
kspace_loc = radial_mask.data
mask = pysap.Image(data=convert_locations_to_mask(kspace_loc, image.shape))

plt.subplot(1, 2, 1)
plt.imshow(np.abs(image), cmap='gray')
plt.title("MRI Data")
plt.subplot(1, 2, 2)
plt.imshow(mask, cmap='gray')
plt.title("K-space Sampling Mask")
plt.show()

$\leadsto$   The code above is to generate the mask, the key part is

radial_mask = get_sample_data("mri-radial-samples")

where we use the package to import the mri-radial-samples already existing.

Generate the kspace

From the 2D brain slice and the acquisition mask, we retrospectively undersample the k-space using a non cartesian acquisition mask. We then grid the kspace to get the gridded solution as a baseline.

# Get the kspace observation values for the kspace locations
fourier_op = NonCartesianFFT(
    samples=kspace_loc,
    shape=image.shape,
    n_coils=cartesian_ref_image.shape[0],
    #implementation='cpu'
    #implementation='finufft' # Use gpuNUFFT for speed
    implementation='gpuNUFFT'
)
'''
NonCartesianFFT : 定义 非笛卡尔傅立叶变换运算符，用于处理 非规则 k-space 采样（如径向、螺旋轨迹等）。
samples=kspace_loc : 指定 k-space 采样点的坐标（即 kspace_loc）。
shape=image.shape : 设置 MRI 图像的形状（匹配输入 image 的尺寸）。
implementation='finufft' : 选择 finufft 作为计算 NUFFT（非均匀傅立叶变换） 的实现方式。
fourier_op 现在是一个 运算符对象，可以用于在非笛卡尔 k-space 进行傅立叶变换。
'''

# kspace_obs = fourier_op.op(cartesian_ref_image)
kspace_obs = fourier_op.op(cartesian_ref_image.data)

Gridded solution

# Gridded solution
grid_space = np.linspace(-0.5, 0.5, num=image.shape[0])
grid2D = np.meshgrid(grid_space, grid_space)
'''
np.linspace(-0.5, 0.5, num=image.shape[0]) : 生成 均匀分布的坐标点，范围在 [-0.5, 0.5] 之间。
image.shape[0] : 设定 坐标点个数，与图像的尺寸保持一致。
np.meshgrid(grid_space, grid_space) : 生成 2D 网格坐标，用于 将非均匀采样数据重新映射到规则网格。
最终结果：grid2D 是 用于插值的 2D 坐标网格
'''

grid_soln = np.asarray([
    gridded_inverse_fourier_transform_nd(kspace_loc, kspace_obs_ch,
                                         tuple(grid2D), 'linear')
    for kspace_obs_ch in kspace_obs
])

'''
gridded_inverse_fourier_transform_nd() : 采用 Gridding 插值方法 进行 k-space 数据的逆傅立叶变换。

kspace_loc : MRI 非笛卡尔 k-space 采样位置（NonCartesianFFT 的 samples）。
kspace_obs : 欠采样 k-space 观测数据（fourier_op.op(image) 计算的结果）。
tuple(grid2D) : 规则网格坐标（转换为 tuple）。
'linear' : 使用 线性插值（Linear Interpolation） 进行 Gridding 重建。
最终结果：grid_soln 是 基于 Gridding 方法的 MRI 影像重建结果。
'''

image_rec0 = pysap.Image(data=np.sqrt(np.sum(np.abs(grid_soln)**2, axis=0)))

plt.imshow(image_rec0, cmap='gray')
base_ssim = ssim(image_rec0, image)
plt.title('Gridded solution : SSIM = ' + str(np.around(base_ssim, 3)))
plt.show()

print('The Base SSIM is : {}'.format(base_ssim))

The Base SSIM is : 0.6173104381198481

Q2  

• Observe and comment the improved image quality on the zero-filled adjoint FFT re construction when using multi-coil data as compared to single coil data.

cartesian_ref_image.data.shape

(32, 512, 512)

cartesian_ref_image.data[0].reshape(1, 512, 512).shape

(1, 512, 512)

cartesian_ref_image.data[0:2, : , :].shape

(2, 512, 512)

fig, axs = plt.subplots(1, 3, figsize=(15, 5))
'''
fourier_op_multi_coil = FFT(mask=mask, shape=image.shape,
                 n_coils=cartesian_ref_image.shape[0])
kspace_obs_multi_coil = fourier_op_multi_coil.op(cartesian_ref_image)
zero_soln_multi_coil = np.linalg.norm(fourier_op.adj_op(kspace_obs), axis=0)
base_ssim_multi_coil = ssim(zero_soln_multi_coil, image)
'''
fourier_op_multi_coil = NonCartesianFFT(
    samples=kspace_loc,
    shape=image.shape,
    n_coils=cartesian_ref_image.shape[0],
    #implementation='cpu'
    #implementation='finufft' # Use gpuNUFFT for speed
    implementation='gpuNUFFT'
)
kspace_obs_multi_coil = fourier_op_multi_coil.op(cartesian_ref_image.data)

# Gridded solution
grid_space = np.linspace(-0.5, 0.5, num=image.shape[0])
grid2D = np.meshgrid(grid_space, grid_space)
grid_soln_multi_coil = np.asarray([
    gridded_inverse_fourier_transform_nd(kspace_loc, kspace_obs_ch,
                                         tuple(grid2D), 'linear')
    for kspace_obs_ch in kspace_obs_multi_coil
])
image_rec0_multi_coil = pysap.Image(data=np.sqrt(np.sum(np.abs(grid_soln_multi_coil)**2, axis=0)))

axs[0].imshow(np.abs(image_rec0_multi_coil), cmap='gray')
base_ssim_multi_coil = ssim(image_rec0_multi_coil, image)
axs[0].set_title('Gridded solution (multi_coil) : SSIM = ' + str(np.around(base_ssim_multi_coil, 3)))
axs[0].axis('off')

'''
fourier_op_single_coil = FFT(mask=mask, shape=image.shape,
                 n_coils=1)
kspace_obs_single_coil = fourier_op_single_coil.op(cartesian_ref_image)
zero_soln_single_coil = fourier_op_single_coil.adj_op(kspace_obs[0])
base_ssim_single_coil = ssim(zero_soln_single_coil, image)
'''

fourier_op_single_coil = NonCartesianFFT(samples=kspace_loc, shape=image.shape, implementation='gpuNUFFT')
kspace_obs_single_coil = fourier_op_single_coil.op(cartesian_ref_image.data[0:1, : , :].reshape(512, 512))

# Gridded solution
grid_space = np.linspace(-0.5, 0.5, num=image.shape[0])
grid2D = np.meshgrid(grid_space, grid_space)
grid_soln_single_coil = gridded_inverse_fourier_transform_nd(kspace_loc, kspace_obs_single_coil,
                                                 tuple(grid2D), 'linear')   # 线性插值

image_rec0_single_coil = grid_soln_single_coil

axs[1].imshow(np.abs(image_rec0_single_coil), cmap='gray')
base_ssim_single_coil = ssim(image_rec0_single_coil, image)
axs[1].set_title('Gridded solution (single_coil) : SSIM = ' + str(np.around(base_ssim_single_coil, 3)))
axs[1].axis('off')

axs[2].imshow(np.abs(image), cmap='gray')
base_ssim_single_coil = ssim(image_rec0_single_coil, image)
axs[2].set_title(f"MRI Data : SSIM = {ssim(image, image)}")
axs[2].axis('off')

plt.tight_layout()
plt.show()

fig, axs = plt.subplots(1, 3, figsize=(15, 5))
'''
fourier_op_multi_coil = FFT(mask=mask, shape=image.shape,
                 n_coils=cartesian_ref_image.shape[0])
kspace_obs_multi_coil = fourier_op_multi_coil.op(cartesian_ref_image)
zero_soln_multi_coil = np.linalg.norm(fourier_op.adj_op(kspace_obs), axis=0)
base_ssim_multi_coil = ssim(zero_soln_multi_coil, image)
'''
fourier_op_multi_coil = NonCartesianFFT(
    samples=kspace_loc,
    shape=image.shape,
    n_coils=cartesian_ref_image.shape[0],
    #implementation='cpu'
    #implementation='finufft' # Use gpuNUFFT for speed
    implementation='gpuNUFFT'
)
kspace_obs_multi_coil = fourier_op_multi_coil.op(cartesian_ref_image.data)

# Gridded solution
grid_space = np.linspace(-0.5, 0.5, num=image.shape[0])
grid2D = np.meshgrid(grid_space, grid_space)
grid_soln_multi_coil = np.asarray([
    gridded_inverse_fourier_transform_nd(kspace_loc, kspace_obs_ch,
                                         tuple(grid2D), 'linear')
    for kspace_obs_ch in kspace_obs_multi_coil
])
image_rec0_multi_coil = pysap.Image(data=np.sqrt(np.sum(np.abs(grid_soln_multi_coil)**2, axis=0)))

axs[0].imshow(np.abs(image_rec0_multi_coil), cmap='gray')
base_ssim_multi_coil = ssim(image_rec0_multi_coil, image)
axs[0].set_title('Gridded solution (multi_coil) : SSIM = ' + str(np.around(base_ssim_multi_coil, 3)))
axs[0].axis('off')

'''
fourier_op_single_coil = FFT(mask=mask, shape=image.shape,
                 n_coils=1)
kspace_obs_single_coil = fourier_op_single_coil.op(cartesian_ref_image)
zero_soln_single_coil = fourier_op_single_coil.adj_op(kspace_obs[0])
base_ssim_single_coil = ssim(zero_soln_single_coil, image)
'''
fourier_op_single_coil = NonCartesianFFT(
    samples=kspace_loc,
    shape=image.shape,
    n_coils=2,
    #implementation='cpu'
    #implementation='finufft' # Use gpuNUFFT for speed
    implementation='gpuNUFFT'
)
kspace_obs_single_coil = fourier_op_single_coil.op(cartesian_ref_image.data[0:2, : , :])

# Gridded solution
grid_space = np.linspace(-0.5, 0.5, num=image.shape[0])
grid2D = np.meshgrid(grid_space, grid_space)
grid_soln_single_coil = np.asarray([
    gridded_inverse_fourier_transform_nd(kspace_loc, kspace_obs_ch,
                                         tuple(grid2D), 'linear')
    for kspace_obs_ch in kspace_obs_single_coil
])
image_rec0_single_coil = pysap.Image(data=np.sqrt(np.sum(np.abs(grid_soln_single_coil)**2, axis=0)))

axs[1].imshow(np.abs(image_rec0_single_coil), cmap='gray')
base_ssim_single_coil = ssim(image_rec0_single_coil, image)
axs[1].set_title('Gridded solution (double_coil) : SSIM = ' + str(np.around(base_ssim_single_coil, 3)))
axs[1].axis('off')

axs[2].imshow(np.abs(image), cmap='gray')
base_ssim_single_coil = ssim(image_rec0_single_coil, image)
axs[2].set_title(f"MRI Data : SSIM = {ssim(image, image)}")
axs[2].axis('off')

plt.tight_layout()
plt.show()

$\leadsto$   When using multi-coil data, the reconstrction image has a sharper contrast compared to single-coil data.

Estimate Sensitivity maps (Smaps)

# Obtain SMaps
Smaps, SOS = get_Smaps(
    k_space=kspace_obs,
    img_shape=fourier_op.shape,
    samples=kspace_loc,
    thresh=(0.01, 0.01),    # The cutoff threshold in each kspace direction
                            # between 0 and kspace_max (0.5)
    min_samples=kspace_loc.min(axis=0),
    max_samples=kspace_loc.max(axis=0),
    mode='gridding',
    method='linear',
    n_cpu=-1,
)

h=3;w=5;
f, axs = plt.subplots(h,w)
for i in range(h):
    for j in range(w):
        axs[i, j].imshow(np.abs(Smaps[3 * i + j]))
        axs[i, j].axis('off')
plt.subplots_adjust(wspace=0,hspace=0)
plt.show()

! pip show pynfft2

Name: pyNFFT2
Version: 1.4.3
Summary: A pythonic wrapper around NFFT
Home-page: https://github.com/ghisvail/pyNFFT.git
Author: Ghislain Vaillant
Author-email: ghisvail@gmail.com
License: 
Location: /home/home/miniconda3/envs/tensor/lib/python3.11/site-packages
Requires: 
Required-by: 

Setup Fourier operators with SENSE

# Setup Fourier Operator with SENSE. This would initialize the
# fourier operators in the GPU.
# For this we need to specify the implementation as gpuNUFFT
# and also pass the Smaps calculated above
#fourier_implementation = 'cpu'
fourier_implementation = 'gpuNUFFT'
#fourier_implementation = 'finufft'

fourier_op_sense = NonCartesianFFT(
    samples=kspace_loc,
    shape=image.shape,
    n_coils=cartesian_ref_image.shape[0],
    smaps=Smaps,
    implementation=fourier_implementation,
)

Q1  

• Use the gpuNUIFFT backend instead of the finufft one to observe a gain in computing speed when estimating the Smaps and computing the SENSE (i.e. zero-order) image.

import time

def  zero_order_reconstruction (implementation_) :
    time_1 = time.time()
    if implementation_ == "gpuNUFFT" :

        fourier_op = NonCartesianFFT(
            samples=kspace_loc,
            shape=image.shape,
            n_coils=cartesian_ref_image.shape[0],
            implementation='gpuNUFFT'
        )

       
    elif implementation_ == "finufft" :

        fourier_op = NonCartesianFFT(
            samples=kspace_loc,
            shape=image.shape,
            n_coils=cartesian_ref_image.shape[0],
            implementation='finufft' # Use gpuNUFFT for speed
        )
       
    # kspace_obs = fourier_op.op(cartesian_ref_image.data)
    kspace_obs = fourier_op.op(cartesian_ref_image.data)
    
    # Gridded solution
    grid_space = np.linspace(-0.5, 0.5, num=image.shape[0])
    grid2D = np.meshgrid(grid_space, grid_space)

    grid_soln = np.asarray([
    gridded_inverse_fourier_transform_nd(kspace_loc, kspace_obs_ch,
                                         tuple(grid2D), 'linear')
    for kspace_obs_ch in kspace_obs
    ])

    image_rec0 = pysap.Image(data=np.sqrt(np.sum(np.abs(grid_soln)**2, axis=0)))
    time_2 = time.time()
    time_zero_order = time_2 - time_1

    base_ssim = ssim(image_rec0, image)

    # Obtain SMaps
    time_3 = time.time()
    Smaps, SOS = get_Smaps(
        k_space=kspace_obs,
        img_shape=fourier_op.shape,
        samples=kspace_loc,
        thresh=(0.01, 0.01),    # The cutoff threshold in each kspace direction
                                # between 0 and kspace_max (0.5)
        min_samples=kspace_loc.min(axis=0),
        max_samples=kspace_loc.max(axis=0),
        mode='gridding',
        method='linear',
        n_cpu=-1,
    )
    time_4  = time.time()
    time_SMaps = time_4 - time_3

    # Setup Fourier Operator with SENSE. This would initialize the
    # fourier operators in the GPU.
    # For this we need to specify the implementation as gpuNUFFT
    # and also pass the Smaps calculated above
    #fourier_implementation = 'cpu'
    time_5  = time.time()
    if implementation_ == "gpuNUFFT" :
        fourier_implementation = 'gpuNUFFT'
    elif implementation_ == "finufft" :
        fourier_implementation = 'finufft'

    fourier_op_sense = NonCartesianFFT(
        samples=kspace_loc,
        shape=image.shape,
        n_coils=cartesian_ref_image.shape[0],
        smaps=Smaps,
        implementation=fourier_implementation,
    )
    time_6  = time.time()
    time_SENSE = time_6 - time_5

    return image_rec0, time_zero_order, base_ssim, time_SMaps, Smaps, time_SENSE

import matplotlib.pyplot as plt
import matplotlib.gridspec as gridspec

for implementation_ in  ["finufft", "gpuNUFFT"] :

    # Perform reconstruction using "finufft" or "gpuNUFFT" 
    image_rec0, time_zero_order, base_ssim, time_SMaps, Smaps, time_SENSE = zero_order_reconstruction(implementation_)

    # Define grid layout: 1 row, 2 columns (left: single image, right: multiple subplots)
    fig = plt.figure(figsize=(12, 6))
    gs = gridspec.GridSpec(1, 2, width_ratios=[1, 2])  # Left takes 1/3, Right takes 2/3

    #   Left: Gridded Solution (Single Image)
    ax1 = plt.subplot(gs[0])
    ax1.imshow(image_rec0, cmap='gray')
    ax1.set_title('Gridded solution : SSIM = ' + str(np.around(base_ssim, 3)) + ' (' + str(implementation_) + ')')
    ax1.axis('off')

    ax1.text(
        0.5, -0.1,  # Position (centered below image)
        f" time_zero_order = {time_zero_order:.3f} sec\n time_SMaps = {time_SMaps:.3f} sec\n time_SENSE = {time_SENSE:.3f} sec",
        ha='center', va='top', fontsize=12, transform=ax1.transAxes
    )

    #   Right: SMaps Subplot (Multiple Images)
    h, w = 3, 5  # Define grid shape for SMaps
    gs2 = gridspec.GridSpecFromSubplotSpec(h, w, subplot_spec=gs[1])  # Nested grid

    for i in range(h):
        for j in range(w):
            ax = plt.subplot(gs2[i, j])
            ax.imshow(np.abs(Smaps[3 * i + j]))
            ax.axis('off')

    # Adjust layout to avoid overlap
    plt.tight_layout()
    plt.show()

/home/home/miniconda3/envs/tensor/lib/python3.11/site-packages/mrinufft/_utils.py:94: UserWarning: Samples will be rescaled to [-pi, pi), assuming they were in [-0.5, 0.5)
  warnings.warn(
[Parallel(n_jobs=-1)]: Using backend LokyBackend with 16 concurrent workers.
[Parallel(n_jobs=-1)]: Done  32 out of  32 | elapsed:    0.3s finished
/home/home/miniconda3/envs/tensor/lib/python3.11/site-packages/mrinufft/_utils.py:94: UserWarning: Samples will be rescaled to [-pi, pi), assuming they were in [-0.5, 0.5)
  warnings.warn(

[Parallel(n_jobs=-1)]: Using backend LokyBackend with 16 concurrent workers.
[Parallel(n_jobs=-1)]: Done  32 out of  32 | elapsed:    0.3s finished
/home/home/miniconda3/envs/tensor/lib/python3.11/site-packages/mrinufft/operators/interfaces/gpunufft.py:146: UserWarning: no pinning provided, pinning existing smaps now.
  warnings.warn("no pinning provided, pinning existing smaps now.")

$\leadsto$   When using gpuNUIFFT ,

• the computation of the zero-order image is slightly faster
• the estimating the Smaps is slightly faster
• and the computation of computing SENSE is significantly slower

compared to the finufft one

FISTA optimization

We now want to refine the zero order solution using a FISTA optimization. The cost function is set to Proximity Cost + Gradient Cost

# Setup the operators
linear_op = WaveletN(wavelet_name='sym8', nb_scale=4)
# regularizer_op = SparseThreshold(Identity(), 4e-7, thresh_type="soft")
regularizer_op = SparseThreshold(Identity(), 1e-8, thresh_type="soft")

# Setup Reconstructor
reconstructor = SelfCalibrationReconstructor(
    fourier_op=fourier_op_sense,
    linear_op=linear_op,
    regularizer_op=regularizer_op,
    gradient_formulation='synthesis',

    lipschitz_cst=fourier_op_sense.impl.get_lipschitz_cst(),  # add

    verbose=1,
)

Lipschitz constant is 16.209672927856445

WARNING: Making input data immutable.

Synthesis formulation: FISTA optimization

We now want to refine the zero order solution using a FISTA optimization. The cost function is set to Proximity Cost + Gradient Cost

x_final, costs, metrics = reconstructor.reconstruct(
    kspace_data=kspace_obs,
    optimization_alg='fista',
    num_iterations=100,
)

image_rec = pysap.Image(data=x_final)
recon_ssim = ssim(image_rec, image)

plt.imshow(np.abs(image_rec), cmap='gray')
plt.title('FISTA Reconstruction\nSSIM = {}'.format(recon_ssim))
plt.show()

 - mu:  1e-08
 - lipschitz constant:  16.209672927856445
 - data:  (512, 512)
 - wavelet:  <mri.operators.linear.wavelet.WaveletN object at 0x7fb2a9afdc90> - 4
 - max iterations:  100
 - image variable shape:  (512, 512)
 - alpha variable shape:  (291721,)
----------------------------------------
Starting optimization...

  0%|          | 0/100 [00:00<?, ?it/s]

100%|██████████| 100/100 [00:14<00:00,  6.84it/s]

 - final iteration number:  100
 - final log10 cost value:  6.0
 - converged:  False
Done.
Execution time:  14.612145459017484  seconds
----------------------------------------

Q3  

• Do the same for the CS-based (synthesis formulation) solutions (FISTA or POGM).

type(cartesian_ref_image)

pysap.base.image.Image

type(cartesian_ref_image.data)

numpy.ndarray

print(kspace_loc.shape)

(32768, 2)

from mri.reconstructors import SingleChannelReconstructor

def reconstruction (coil):


    if coil == "multi_coil" :

        fourier_op_sense = NonCartesianFFT(
            samples=kspace_loc,
            shape=image.shape,
            n_coils=cartesian_ref_image.shape[0],
            smaps=Smaps,
            implementation='gpuNUFFT',
        )

        fourier_op = NonCartesianFFT(
            samples=kspace_loc,
            shape=image.shape,
            n_coils=cartesian_ref_image.shape[0],
            implementation='gpuNUFFT'
        )

        kspace_obs = fourier_op.op(cartesian_ref_image.data)

        # Setup the operators
        linear_op = WaveletN(wavelet_name="sym8", nb_scales=4)
        regularizer_op = SparseThreshold(Identity(), 1e-8, thresh_type="soft")

        # Setup Reconstructor
        reconstructor = SelfCalibrationReconstructor(
            fourier_op=fourier_op_sense,
            linear_op=linear_op,
            regularizer_op=regularizer_op,
            gradient_formulation='synthesis',

            lipschitz_cst=fourier_op_sense.impl.get_lipschitz_cst(),  # add

            verbose=1,
        )

        x_final, costs, metrics = reconstructor.reconstruct(
            kspace_data=kspace_obs,
            optimization_alg='fista',
            num_iterations=100,
        )

        image_rec = pysap.Image(data=x_final)
        recon_ssim = ssim(image_rec, image)

    elif coil == "single_coil" :

        fourier_op = NonCartesianFFT(samples=kspace_loc, shape=image.shape, implementation='gpuNUFFT')
        
        kspace_obs = fourier_op.op(cartesian_ref_image.data[0])

        # Setup the operators
        linear_op = WaveletN(wavelet_name='sym8', nb_scale=4)
        regularizer_op = SparseThreshold(Identity(), 1e-8, thresh_type="soft")

        # Setup Reconstructor
        reconstructor = SingleChannelReconstructor(
            fourier_op=fourier_op,                                              # the **Non Cartesian  Fourier transform** (NonCartesianFFT) operator only for the sampled locations ( `kspace_loc` )
            linear_op=linear_op,                                                # Symmlet-8 wavelet transform with 4 scales for multi-scale sparse representation.
            regularizer_op=regularizer_op,                                      # Identity() < - > identity linear transformation, \lambda = 6 * 1e-7, Soft Thresholding < - > L1 norm
            gradient_formulation='synthesis',                                   # using Synthesis Formulation
            verbose=1,
        )

        image_rec, costs, metrics = reconstructor.reconstruct(
            kspace_data=kspace_obs,                                              # kspace_obs = fourier_op.op(image)
            optimization_alg='fista',                                            # fista
            num_iterations=200,
        )

        recon_ssim = ssim(image_rec, image)

    return image_rec, recon_ssim

fig, axs = plt.subplots(1, 3, figsize=(15, 5))


FISTA_soln_multi_coil, FISTA_ssim_multi_coil = reconstruction("multi_coil")
axs[0].imshow(np.abs(FISTA_soln_multi_coil), cmap='gray')
axs[0].set_title('Iterative FISTA Reconstruction (multi_coil) : SSIM = ' + str(np.around(FISTA_ssim_multi_coil, 3)))
axs[0].axis('off')

FISTA_soln_single_coil, FISTA_ssim_single_coil = reconstruction("single_coil")

axs[1].imshow(np.abs(FISTA_soln_single_coil), cmap='gray')
axs[1].set_title('Iterative FISTA Reconstruction (single_coil) : SSIM = ' + str(np.around(FISTA_ssim_single_coil, 3)))
axs[1].axis('off')

axs[2].imshow(np.abs(image), cmap='gray')
base_ssim_single_coil = ssim(image_rec0_single_coil, image)
axs[2].set_title(f"MRI Data : SSIM = {ssim(image, image)}")
axs[2].axis('off')

plt.tight_layout()
plt.show()

WARNING: Making input data immutable.
WARNING: Making input data immutable.

Lipschitz constant is 16.209671020507812
 - mu:  1e-08
 - lipschitz constant:  16.209671020507812
 - data:  (512, 512)
 - wavelet:  <mri.operators.linear.wavelet.WaveletN object at 0x7fb2f7cc2290> - 4
 - max iterations:  100
 - image variable shape:  (512, 512)
 - alpha variable shape:  (291721,)
----------------------------------------
Starting optimization...

100%|██████████| 100/100 [00:14<00:00,  6.97it/s]
WARNING: Making input data immutable.

 - final iteration number:  100
 - final log10 cost value:  6.0
 - converged:  False
Done.
Execution time:  14.347829079022631  seconds
----------------------------------------

WARNING: Making input data immutable.

Lipschitz constant is 17.830640220642092
 - mu:  1e-08
 - lipschitz constant:  17.830640220642092
 - data:  (512, 512)
 - wavelet:  <mri.operators.linear.wavelet.WaveletN object at 0x7fb13c585ed0> - 4
 - max iterations:  200
 - image variable shape:  (512, 512)
 - alpha variable shape:  (291721,)
----------------------------------------
Starting optimization...

100%|██████████| 200/200 [00:10<00:00, 19.41it/s]

 - final iteration number:  200
 - final log10 cost value:  6.0
 - converged:  False
Done.
Execution time:  10.306300255993847  seconds
----------------------------------------

$\leadsto$   When using multi-coil data, the quality of iterative FISTA reconstruction is better to that of single-coil data.

POGM reconstruction

x_final, costs, metrics = reconstructor.reconstruct(
    kspace_data=kspace_obs,
    optimization_alg='pogm',
    num_iterations=200,
)

image_rec = pysap.Image(data=np.abs(x_final))
recon_ssim = ssim(image_rec, image)

plt.imshow(np.abs(image_rec), cmap='gray')
recon_ssim = ssim(image_rec, image)
plt.title('POGM Reconstruction\nSSIM = {}'.format(recon_ssim))
plt.show()

 - mu:  1e-08
 - lipschitz constant:  16.20966911315918
 - data:  (512, 512)
 - wavelet:  <mri.operators.linear.wavelet.WaveletN object at 0x7fb2f7a0d2d0> - 4
 - max iterations:  200
 - image variable shape:  (1, 512, 512)
----------------------------------------
Starting optimization...

100%|██████████| 200/200 [00:28<00:00,  7.04it/s]

 - final iteration number:  200
 - final log10 cost value:  6.0
 - converged:  False
Done.
Execution time:  28.399563284998294  seconds
----------------------------------------

Analysis formulation: Condat-Vu reconstruction

#linear_op = WaveletN(wavelet_name="sym8", nb_scales=4)
linear_op = WaveletUD2(
    wavelet_id=24,
    nb_scale=4,
)
# regularizer_op = SparseThreshold(Identity(), 4e-7, thresh_type="soft")
regularizer_op = SparseThreshold(Identity(), 1e-9, thresh_type="soft")

# Setup Reconstructor
reconstructor = SelfCalibrationReconstructor(
    fourier_op=fourier_op_sense,
    linear_op=linear_op,
    regularizer_op=regularizer_op,
    gradient_formulation='analysis',
    verbose=1,
)

WARNING: Making input data immutable.

Lipschitz constant is 11.979288864135743

x_final, costs, metrics = reconstructor.reconstruct(
    kspace_data=kspace_obs,
    optimization_alg='condatvu',
    num_iterations=200,
)
image_rec = pysap.Image(data=np.abs(x_final))
plt.imshow(np.abs(image_rec), cmap='gray')
recon_ssim = ssim(image_rec, image)
plt.title('Condat-Vu Reconstruction\nSSIM = {}'.format(recon_ssim))
plt.show()

WARNING: <class 'mri.operators.linear.wavelet.WaveletUD2'> does not inherit an operator parent.

 - mu:  1e-09
 - lipschitz constant:  11.979288864135743
 - tau:  0.15369599612293275
 - sigma:  0.5
 - rho:  1.0
 - std:  None
 - 1/tau - sigma||L||^2 >= beta/2:  True
 - data:  (512, 512)
 - wavelet:  <mri.operators.linear.wavelet.WaveletUD2 object at 0x7fb2f7bee610> - 4
 - max iterations:  200
 - number of reweights:  0
 - primal variable shape:  (512, 512)
 - dual variable shape:  (2621440,)
----------------------------------------
Starting optimization...

100%|██████████| 200/200 [04:16<00:00,  1.28s/it]

 - final iteration number:  200
 - final cost value:  1000000.0
 - converged:  False
Done.
Execution time:  258.0937590090034  seconds
----------------------------------------